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We investigate the Anderson transition found in the spectrum of the Dirac operator of Quantum 
Chromodynamics (QCD) at high temperature, studying the properties of the critical quark eigen¬ 
functions. Applying multifractal finite-size scaling we determine the critical point and the critical 
exponent of the transition, finding agreement with previous results, and with available results for 
the unitary Anderson model. We estimate several multifractal exponents, finding also in this case 
agreement with a recent determination for the unitary Anderson model. Our results confirm the 
presence of a true Anderson localization-delocalization transition in the spectrum of the quark Dirac 
operator at high-temperature, and further support that it belongs to the 3D unitary Anderson model 
class. 
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I. INTRODUCTION 

The Anderson metal-insulator transition is a genuine 
quantum phase transition, which has been widely investi¬ 
gated in condensed matter physics since the seminal pa¬ 
per of AndersorP. In the past years Anderson transitions 
were found in a wide range of physical s yst ems, such as 
ultrasound in disordered elastic networks^, light in dis¬ 
ordered photonic lattices in the transverse directiorP, or 
in an ultracold atomic system in a disordered laser trapP. 

A characteristic feature of Anderson transitions is the 
rich multifractal structure of critical eigenstates, which 
has been the subject of intense research activity in re¬ 
cent years (see Ref. [6] for a review). Direct signs of 
multifractals at the metal-insulator transition point have 
been observed experimentally in dilute magnetic semi¬ 
conductor^. Multifractality can moreover influence the 
behavior of various systems near criticality in different 
ways. For example, the large overlap of multifractal 
wave-functions can increase the superconducting critical 
temperatur^. The multifractality of the local density of 
states may induce a new phase because of the presence 
of local Kondo effects induced by local pseudogaps at the 
Fermi energjP. 

The simplest model displaying a metal-insulator tran¬ 
sition is the Anderson model, which describes non¬ 
interacting fermions in a disordered crystal. Disorder is 
usually introduced through a random on-site potential, 
while hopping elements are fixed (up to a random phase 
or SU(2) rotation). In this case the system belongs to one 
of the Wigner-Dyson (WD) symmetry classes depending 
on the global symmetries of the system. On the other 
hand, systems with vanishing on-site terms and random 
hopping terms, if the lattice is bipartite, possess an ad¬ 
ditional chiral symmetry and belong to one of the chiral 


WD classed. 

Quite surprisingly, an Anderson transition has been 
shown to take place also in the spectrum of the Dirac 
operator in Quantum Ghromodynamics (QGD) at high 
temperatur^iSHIll (see Ref. [TS] for a review). QCD is 
the quantum field theory governing strong interactions 
at the microscopic level, and operates on length and en¬ 
ergy scales vastly different from the ones usually encoun¬ 
tered in condensed matter physics. QCD is a non Abelian 
gauge theory, describing the interactions of quarks, which 
are fermions, and gluons, which are the vector bosons 
of the SU{3) gauge symmetry. Although these are the 
fundamental degrees of freedom, they do not appear in 
the spectrum of the theory at low temperatures, which 
contains only hadrons, i.e., bound states of quarks and 
gluons, due to the phenomenon of confinement. However, 
at a (pseudo)critical temperature, Tc, strongly interact¬ 
ing matter undergoes a crossover to the so-called quark- 
gluon-plasma phase, and at high temperatures quarks 
and gluons are deconfined. This transition is accompa¬ 
nied by the restoration of the approximate chiral symme¬ 
try that is spontaneously broken at low temperature^^. 

Contributions of quarks to observables, as well as all 
quark correlation functions, are entirely encoded in the 
eigenvalues and the eigenvectors of the Dirac operator 
in the background of a non Abelian gauge field. Phys¬ 
ical quantities are then obtained after averaging over 
the gauge field configurations with the appropriate path- 
integral measure. In this respect, the eigenmodes of the 
Dirac operator can be formally treated as the eigenstates 
of a random “Hamiltonian”, with disorder provided by 
the gauge field fluctuations. Among the eigenmodes, a 
prominent role is played by the low-lying ones: for ex¬ 
ample, they give the most important contributions to 
the quark correlation functions, and determine the fate 
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FIG. 1. Eigenvectors of the Dirac operator in lattice QCD at T « 2.6Tc (a) in the insulating regime, (b) at criticality, and (c) 
in the metallic regime. Dot sizes are proportional to with proportionality factors tuned independently for each subfigure 

to improve visualization. The spatial system size is L = 56 (in lattice units) for all subfigures. Coloring is determined by the 
value of the x coordinate. 


of chiral symmetry through the Banks-Casher relatiorfi^. 
The low end of the spectrum looks completely different 
in the hadronic and in the quark-gluon-plasma phase. 
At low temperatures, the density of states is finite near 
the origin, and both low-lying and bulk eigenmodes are 
extended throughout the system. In contrast, at high 
temperatures, above T^, the density of states vanishes at 
the origin, and the low-lying eigenmodes are localized on 
the scale of the inverse temperature, while higher up in 
the spectrum, beyond a temperature-dependent c ritical 
“energy”, EdT), the eigenmodes are again extendecP^EH. 
The temperature dependence of the “mobility edge” was 
investigated in Ref. [HI in which it was found that Ec{T) 
extrapolates to zero at a temperature compatible with T^. 
Typical Dirac eigenmodes in the localized, critical and 
delocalized regimes are shown in Fig. The transition 
in the spectrum from localized to delocalized eigenmodes 
has been shown to be a second-order phase transition, 
with critical exponent compatible with the one found in 
the three-dimensional unitary Anderson mode^. 

It is rather surprising at first that the Anderson transi¬ 
tion in the high-temperature QCD Dirac spectrum seems 
to belong to the same universality class as that of the 
three-dimensional unitary Anderson model. From the 
point of view of statistical systems, QCD at a finite tem¬ 
perature T is in fact a four-dimensional Euclidean model, 
with the “temporal” dimension compactified on a circle 
of length 1/T. However, it has been argued that high- 
temperature QCD is an effectively three-dimensional dis¬ 
ordered system with on-site disorder, the strength of 
which is set by the temperatur (j^^l^°l While this makes it 
more plausible that the two models actually belong to the 
same universality class, it does not make it less impor¬ 
tant to look for further evidence. In this respect, finding 


the same multifractal structure in the critical eigenstates 
would give strong support to the claim of Ref. [HI and 
so to the broader universality of the critical properties 
of Anderson transitions. The study of this multifractal 
structure is precisely the aim of this work. 

The plan of the paper is the following. In section [IT] 
we give a brief discussion of multifractality, and of the 
method of multifractal finite-size scaling (MFSS). In sec¬ 
tion |III| we describe in some detail the Dirac operator 
and the numerical simulations of QCD employed in this 
paper. In section IV we study the correlations between 
eigenvectors of the Dirac operator around the critical en¬ 
ergy, both for comparison to the 3D unitary Anderson 
model, and for their appropriate treatment in the sta¬ 
tistical analysis. In section [V] we discuss the results of 
MFSS for the eigenvectors of the Dirac operator. Finally, 
in section m we state our conclusions. 


II. FINITE-SIZE SCALING LAWS FOR 
GENERALIZED MULTIFRACTAL EXPONENTS 

In this section we briefly review wave-function mul¬ 
tifractality and the technique of multifractal finite-size 
scaling. The wave-function of ^ particle in is 

naturally associated to a local probability distribution, 
namely giving the probability to find the par¬ 

ticle in an infinitesimal neighborhood of x. For smooth 
wave functions, the probability to find the particle in 
a small but finite neighborhood of x of size r scales as 
~ r'^. For fractal wave functions, this probability scales 
as ~ r“, where a < c? is called the fractal dimension. For 
strongly fluctuating wave-functions, however, this prob¬ 
ability scales in general as ~ with an x-dependent 
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power a(x) called the local dimension. In turn, points 
with the same local dimension, a(x) = a, constitute 
a subset of characterized by its own fractal dimen¬ 
sion, which generally depends on a. The wave-function 
therefore defines not one, but many different fractals, 
and is therefore said to be multifractal. Multifractal 
wave-functions are strongly fluctuating on every length- 
scale, and their characterization requires an infinite num¬ 
ber of fractal dimensions, called multifractal exponents 
(MFEs). An example of a multifractal wave-function is 
shown in Fig. [^b). 

Multifractality is a known feature of critical eigenfunc¬ 
tions at the Anderson metal-insulator transitiorP, that 
can be studied by means of multifractal finite-size scal¬ 
ing (MFSSJPl. In recent high-precision calculation^^HMl^ 
MFSS has been successfully employed to determine the 
MFEs of critical eigenfunctions, as well as to obtain a 
more precise estimate of the critical disorder and of the 
critical exponents, for Anderson models in different sym¬ 
metry classes. In this work we want to perform a sim¬ 
ilar MFSS analysis to study the Anderson localization- 
delocalization transition in the spectrum of the Dirac op¬ 
erator in QCD. 

In the remainder of this section we describe MFSS in 
some detail. Our methods and notations are essentially 
the same as in Ref. [231 to which we refer the reader for 
a more detailed discussion. There is however one impor¬ 
tant difference, concerning the way in which the transi¬ 
tion is approached. In Ref. 1231 the transition was studied 
by looking at wave functions at the band center and vary¬ 
ing the amount of disorder, W. In QCD the amount of 
disorder is effectively set by the temperature, and it is 
more convenient to keep it fixed and study the transition 
as a function of energy, E, by looking at wave-functions 
near the mobility edge, Ec- Therefore, W has been re¬ 
placed by E in the expressions of Ref. [231 

Let us consider a d-dimensional cubic lattice of linear 
size L, and a critical eigenfunction of a random Hamilto¬ 
nian, ijj{x), defined on the lattice sites x and normalized 
to 1. We can divide the lattice into smaller boxes of lin¬ 
ear size £, and compute the probability corresponding to 
the /c-th box as 


Mfe = XI ( 1 ) 

xGboxfc 


where the sum runs over the lattice sites contained in the 
fc-th box. The generalized inverse participation ratios 
(GIPRs) are the moments of the box probability. The 
GIPRs and their derivatives read 


^9 = X 


k^l 




dRq 

dq 


x-‘‘ 


( 2 ) 


where X = and the sum runs over all the boxes of 
size £. For small A, the averages of Rq and Sq over disor¬ 
der realizations follow a power-law behavior as a function 


of A, which leads one to define the following exponents: 


D, 


lim - 

A-j-o q — 


In (Rq) 

1 InA 


Xq = lim — 
^ A^o [R, 


{S, 


' In A 


( 3 ) 


Dq and Uq are generalized fractal dimensions, usually 
referred to as multifractal exponents (MFEs). One can 
similarly define MEEs for localized and delocalized states 
by substituting critical eigenfunctions with localized or 
delocalized eigenfunctions in Eq. 0- In the delocal¬ 
ized/metallic part of the spectrum, states extend over the 
whole lattice, so their effective size grows proportionally 
to the volume, thus leading to D™®* = d. On the other 
hand, in the localized/insulating regime, states are ex¬ 
ponentially localized, so that their effective size does not 
change with the system size, resulting in = 0 for 

q > 0, and D™'* = oo for g < 0. At criticality, E = Ac, 
the eigenstates are instead expected to be multifractal, 
with nontrivial, g-dependent Dq and aq. 

This jump of the MEEs at the critical point happens 
only in an infinite system. The main idea of MFSS is to 
use the MFEs as order parameters for finite size-scaling. 
In order to do that we have to define the finite size version 
of the MEEs at a given energy. 


D^q-%E,LJ) = 


jS,) 

{Rg)lnX^ 

I In {Rq) 
q — 1 In A 


( 4 ) 

( 5 ) 


where it is understood that wave-functions of energy 
around E are used on the right-hand side, and where 
the superscript ens is to remind the reader that one has 
to perform ensemble averaging over the different disorder 
realizations, ctq and Dq are called generalized multifrac¬ 
tal exponents (GMFEs). Every GMFE approaches the 
value of the corresponding MEE at the critical point, 
E = Ec, only in the limit A —>■ 0. One can also define 
typical MFEs, 


alyp{E,L,i) 

blyp{E,L,£) 


\Rq/ InA’ 
1 (InRg) 
q — 1 In A 


( 6 ) 

( 7 ) 


which can be used as well in a finite-size scaling analysis. 
However, as we said above, MFEs are defined through 
ensemble averaging in principle [see Eq. ([^], and when 
computing MFEs in Sec. |V] we use ensemble averaged 
quantities only. 

In the renormalization group language, the Anderson 
transition is characterized by a single relevant operatoi!^, 
and so in the vicinity of the critical point one can derive 
scaling laws for the GMEEs, which can be summarized 
in a single equation, using a common letter, G, for the 
GMEEs: 


Gq{E,L,i)=Gq+:^gq{^,"- 


L l_ 

VI 


( 8 ) 
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At the critical point, the localization length diverges as 
^ - [q{E - Ec)]-"', where giE - E^) - E - E^ for {E - 
Ec)/Ei. 1. The system sizes employed in this paper, 
however, are not big enough to justify the use of one- 
parameter scaling, and so we included the contribution 
of an irrelevant operator, 77 = g^E — Ec), which leads us 
to write 


Gq{E,L,l) — Gq -\- 


1 

In A - 
gt-yg" 


[qL- 


(9) 


A second irrelevant term, proportional to L~y , is ex¬ 
pected to be less important and will be neglected in the 

analysi^2lI13]_ 

Fits to the numerical data are performed by expanding 
g^ and g" in the variables qL^ and gi^ up to order 
and riir, respectively. The number of parameters there¬ 
fore grows as ^ -I- Moreover, g and 77 must also 

be expanded in powers of E — Ec up to order Ug and tt,^, 
which further increases the number of fitting parameters. 
The fit provides all the physically interesting quantities, 
namely the critical point, E^ the critical exponent, 
the irrelevant exponent 77 , and the MFE, Gq. 

A simpler fit can be performed by setting X = ^/L to 
a fixed value. In this case, dropping A from the notation, 
we can write 


Gq{E, L) = gq (^1) = gi (pL^) + gL-yg" , 

( 10 ) 

having absorbed Gq and the factor A *'/lnA into new 
functions . The main advantage is that since g'^J'’’^ 
are now single-variable functions, the number of expan¬ 
sion parameters grows only as ~ + nir- On the other 

hand, with this method one can determine only Ac, iz, 
and 7/, while the value of the MFE, G^, cannot be ob¬ 
tained. 


III. PROPERTIES OF THE DIRAC OPERATOR 
AND DETAILS OF THE SIMULATIONS 

In this section we give the relevant details about the 
Dirac operator and QCD, and about how the QCD Dirac 
spectrum can be studied by means of numerical simula¬ 
tions. The continuum Euclidean Dirac operator is 

4 

-D(A) = ^7^(9^(11) 

where 7 ^ are the Euclidean Dirac matrices, g is the cou¬ 
pling constant, and Ag is the non Abelian gauge field. 
More precisely, A)(i“ is a Hermitian 3x3 ma¬ 
trix, where A“ = A“(a;) = is real and the sum 

runs over the generators of SU (3). The Dirac operator 
is thus anti-Hermitian, so admitting a straightforward in¬ 
terpretation as (z times) the Hamiltonian of a quantum 


system. The Dirac operator is a chiral operator with the 
following structure in spinor space, 

° = ")> 

with W a complex matrix with no further symmetrji^. 
As a random matrix model, the Dirac operator in a ran¬ 
dom gauge field belongs therefore to the chiral unitary 
class. Chiral symmetry is expressed by the anticommuta¬ 
tion relation { 75 , D} = 0 , which implies that the nonzero 
eigenvalues come in pairs EiEn- It is thus sufficient to 
consider the positive part of the spectrum only. 

The partition function of QCD at temperature T can 
be expressed as a functional integral, 

^QCD = AdA] J^det[iA(A)-b m/], (13) 


with the constraint Ag{x,l/T) = A^(a:, 0). The prod¬ 
uct is over the six different types of quarks (“flavors”), 
with ruf the mass of quark /. Here S'g[A] is a posi¬ 
tive functional of the gauge field, which together with 
the determinants provides the probability distribution of 
the disorder, i.e., of the gauge field configurations. Nu¬ 
merical simulations of QCD require the discretization of 
Eq. (13) on a finite lattice. Eor a review of lattice QCD 
see, e.g.. Ref. nzi While the discretization of the gauge 
fields poses no particular problem, and can be performed 
preserving exact gauge invarianc^^, fermion fields are 
known to be more problematic, and the discretization of 
the Dirac operator spoils some of the properties of its con¬ 
tinuum counterpart. Nevertheless, the discretization that 
we employed, namely staggered fermion^^, preserves the 
anti-Hermiticity and the symmetry of the spectrum with 
respect to the origin, and moreover preserves the chiral 
unitary symmetry clas^^. 

It must be noted at this point that in the case of the 
Anderson model, chiral and non-chiral symmetry classes 
differ only in their properties near the band centeJ^, i.e., 
A = 0 , while the properties of the bulk of the spectrum 
are similar. Eor example, the authors of Ref. 1301 found 
Wigner-Dyson statistics in the bulk spectrum of a three- 
dimensional chiral orthogonal disordered model. More¬ 
over, even the critical exponent of the orthogonal and of 
the chiral orthogonal class turn out to be the same, up 
to very high numerical precisiorP^. We expect the same 
to be true for the multifractal exponents. 

Let us now describe the numerical setting in some de¬ 
tail. QCD is discretized on a periodic hypercubic lat¬ 
tice Xg S Z, of spatial extent L in each direction and 
temporal extent Lt. The gauge fields A^ are replaced 
by corresponding gauge links, i.e., parallel transporters 
along each link of the lattice, which are elements of the 
gauge group, SU{3). The functional Sg is discretized 
and expressed in terms of the gauge links, and the inte¬ 
gration over gauge fields is replaced by the integration 
with the Haar measure over gauge links, i.e., over the 
gauge-group valued variables on the links. Finally, the 
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continuum Dirac operator is replaced by the staggered 
Dirac operator, which reads 

1 4 

= 2 ^ [4+A.yt^/.(a;) - 5a,-ji,yUl{x - A)] , 

(14) 

with 7 y^(a;) = (—, and Uy{x) G SU{S) the gauge 
link connecting the lattice site x to the neighboring site 
along direction fi. The staggered Dirac operator carries 
only spacetime and color indices, i.e., it has no spinorial 
structure. The eigenvalue equation must 

be supplemented with the antiperiodic boundary condi¬ 
tion x(x,Lt) = —x(x, 0 ) for the quark eigenfunction. 

As we have already remarked, the Dirac operator can 
be viewed as a random Hamiltonian, with disorder pro¬ 
vided by the fluctuations of the gauge fields, and dis¬ 
tributed according to the Boltzmann weight appearing 
in the partition function. In its discretized version, the 
Dirac operator is a large sparse matrix, with nonzero ran¬ 
dom elements only in the off-diagonal, nearest-neighbor 
hopping terms, which depend on the parallel transporter 
on the corresponding link of the lattice. This resembles 
an Anderson model with off-diagonal disorder, although 
here the fluctuations of the gauge links are correlated, 
rather than independent. However, since the theory has 
a mass gap, correlations decrease exponentially with the 
distance. Moreover, the strong correlation between the 
different time-slices makes the model effectively three- 
dimensional, with the fluctuations of the temporal links 
ac ting e ffectively as a three-dimensional diagonal disor- 
de j — I—I. The size of the gauge field fluctuations are de¬ 
termined by the temperature, which therefore is expected 
to play the same role as the amount of disorder in the 
Anderson model. This is confirmed by the fact that the 
temperature governs the position of the mobility edge. 

In the present work we have studied the spectrum of 
the Dirac operator by generating gauge link configura¬ 
tions, i.e., realizations of disorder, by means of Monte- 
Carlo methods. Numerical calculations were done on a 
GPU cluster. In our simulations we have included only 
the three lightest flavors (up, down, and strange), with 
equal masses for the up and down quark. For many pur¬ 
poses, this is a good approximation of the real world. The 
lattice spacing in physical units was set to a = 0.125 fm 
and the temporal size was fixed to L* = 4, resulting in 
the temperature T « 2.6 Tc, well above the crossover 
temperature (see Refs. [TTVfT^ for more details). Techni¬ 
cal details about the numerical implementation and the 
scale-setting procedure can be found in Refs. 15^ and 1^ 
We have computed the eigenpairs of the Dirac operator 
from the smallest eigenvalue up to the upper end of the 
critical region, on lattices of spatial sizes in the range 
L = 24 — 56 (in lattice units). A detailed list is reported 
in Tab. |T] along with the corresponding number of sam¬ 
ples. 

The three-dimensional box probability, Eq. 0 , re¬ 
quired for the multifractal analysis, was constructed 
as follows. To have a gauge-invariant description we 


system size (L) 

number of samples 

24 

41517 

28 

20548 

32 

19250 

36 

14869 

40 

8812 

44 

5242 

48 

7008 

56 

3107 


TABLE I. System sizes and corresponding number of gauge 
configurations used in this work. 

(a) (b) 




FIG. 2. Correlations, Eq. ( |17[ ), between (a) critical eigen¬ 
functions of the unitary Anderson model for W — 18.37 and 
system size L = 10, and (b) eigenfunctions of the QCD Dirac 
operator at T ~ 2.6 Tc in the critical regime, 0.32 < E < 0.35, 
for different system sizes. 


summed over the color components, labelled by c. More¬ 
over, due to the strong correlation between the lattice 
time-slices, the eigenvectors of the Dirac operator look 
qualitatively the same on each of them, so we can also 
sum over the time-slices, t. The squared amplitude 
|'!/'(^)P is then defined as |'!/'(^)P = c IXc(u 0 p 5 ^^id 
provides the basic three-dimensional spatial probability 
distribution, from which the box probability distribution 
is then obtained in the usual way. 


IV. CORRELATIONS BETWEEN 
EIGENVECTORS 

In this section we investigate the correlations between 
different eigenvectors of the Dirac operator in a given 
gauge configuration. Our motivation is twofold. On the 
one hand, we want to compare the eigenvector correla¬ 
tions in QCD with the ones found in the unitary Ander¬ 
son model. On the other hand, these correlations have to 
be properly taken into account when fitting the numer¬ 
ical data to determine the various critical quantities, as 
we do in Sec. El 

Cuevas and Kravtso\Ell showed that in the Ander¬ 
son model there are non-negligible correlations between 
eigenvectors. Similar correlations are therefore expected 
also in other disordered systems, like the one under 
consideration. The relevant quantities are the density- 
density correlations, which are defined in terms of the 
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overlap integral, which for the i-th and j-th eigenfunc¬ 
tions reads 

^2 = J ■ ( 15 ) 


In the case of QCD, has the meaning discussed 

above at the end of Sec. |III| One then defines the joint 
probability distribution of ol Ih® energy differ¬ 

ence between eigenstates, 

P(a;, k) = S{Ei - Ej - uj)5{KI^ - . (16) 

To characterize the average behavior of the overlap inte¬ 
gral as a function of energy, its conditional expectation 
value is the natural choice. 


^ _ JdkkP{uj,k) 

^ ’ JdkP{uj,k) • 


(17) 


The quantity 0(0;) is expected to be of order 1/N along 
the whole spectrum, where N = is the volume of the 
system. Indeed, for two delocalized states « 1/lV, 
while for two localized states K 2 is nonzero only if they 
happen to be in the same region, in which case it is of 
order 1, and the probability that this happens is of order 
1/N. 

Fig.j^a) shows the eigenvector correlation C{uj) in the 
unitary Anderson model at criticality. One can see a large 
enhancement of the correlation at small w, and decreas¬ 
ing behavior with growing energy separations, which is 
similar to the results of Ref. [34] for the orthogonal Ander¬ 
son model. Examining the same correlation for critical 
eigenfunctions in QCD, we also find an enhancement at 
small energy separations, see Fig. [^b). In the critical 
regime the behavior of the two systems is very similar, 
and even the approximate exponent of the power-law de¬ 
cay is close to 0.5 in both cases. This is a nice example of 
the similarity of the two models, and in Sec.j^we present 
further similarities in more detail. 


V. MFSS FOR THE EIGENVECTORS OF THE 
DIRAC OPERATOR 

In this section we would like to characterize the Ander¬ 
son phase transition in the spectrum of the Dirac opera¬ 
tor of QCD in the frame of the MFSS, described in Sec.jnj 
As discussed at the end of Sec. m a three-dimensional 
spatial probability distribution was calculated from the 
eigenvectors. From that, the GMFEs and Dg were 
then computed according to Eqs. Q-Q. More pre¬ 
cisely, we chose 26 values of energy, Ei, in the range 
E e [0.32,0.35], and for the Tth energy value and the 
A:-th gauge configuration we computed and ac¬ 
cording to Eq. In order to decrease the numerical 
noise we averaged over all the eigenvectors in an energy 
range of width AE = 0.0012 around Ei. The GMFEs 


aq{Ei, L, tj and Dq{Ei, L, i) are then obtained by averag¬ 
ing and Sqi over the index k, i.e., over configurations, 
or in other words, over different realizations of disorder. 

An example of the resulting GMFEs at fixed A = 0.125 
is depicted in Fig. As the system size grows, the curves 
shift to opposite directions on the two sides of the tran¬ 
sition. At low energy they shift down, indicating a local¬ 
ized phase, while at high energy they shift up, suggest¬ 
ing a metallic phase, as expected. In between, the curves 
should cross at a common point, corresponding to the 
critical energy, but due to finite size effects originating 
from the irrelevant terms this is true only approximately. 

Data were then fitted with the scaling laws Eqs. 
and ( [l0| ), minimizing the quantity y^/(A’d/ —1), using the 
MINOTT librarjf^. Here N^f is the number of degrees of 
freedom and is the distance between the numerical 
data, Ui, and the fitting function, fi, in the appropriate 
metric, i.e., 




where C is the covariance matrix of the data points. In 
the light of the results of Sec. |IV[ which show that there 
are strong correlations between eigenvectors in a given 
gauge configuration, strong correlations are also expected 
among GMFEs at different energies, and so the inclusion 
of correlations in the fitting procedure is necessary to 
obtain accurate results. The error bars of the best fit 
parameters were estimated by Monte-Garlo simulation, 
generating Nmc = 100 sets of synthetic data, distributed 
according to a Gaussian distribution with means equal to 
the raw data points and covariance matrix equal to the 
covariance matrix of the sample. We then determined 
the error bars from the distributions of the resulting fit 
parameters, choosing the 95% confidence level. 

In order to perform best fits, the scaling laws Eqs. © 
and (10) need to be expanded in powers of if — E^, and 
this requires to set the expansion orders n^/ir of the rel¬ 


evant/irrelevant scaling term Q/J"., as well as the expan¬ 
sion orders rig and of g and g. Since the relevant 
operator is more important than the irrelevant one we 
always used n.^ > ?T-ir and rig > rig. We then repeated 
the fit for several choices of the expansion orders. 

The quality of the best fits was judged according to 
two criteria. The first criterion was how close the ra¬ 
tio X^ /{Ndf — 1) approached unity, and only fits with 
X^/{Ndf — 1) « 1 were considered acceptable. The sec¬ 
ond criterion was stability against changing the expan¬ 
sion orders, in order to keep under control the system¬ 
atic effects due to the truncation of the scaling function. 
We estimated the systematic error due to truncation as 
twice the standard deviation of the critical parameters, in 
the sample comprising the stable fits and the essentially 
equivalent ones obtained by increasing or lowering the ex¬ 
pansion orders. The factor of 2 is required by consistency 
with the 95% confidence level chosen for the statistical 
error. 
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(a) (b) 




FIG. 3. GMFEs, (a) Dq"i and (b) a{^Q, at fixed A = 0.125. Dots are the raw data, and the solid red line is the best fit obtained 
by MFSS. Insets show the scaling functions on a log-log scale, after subtracting the irrelevant term. Error bars are not shown 
in the insets for visual clarity. 


We first performed the MFSS at fixed A, as described 
in Sec. [ITj both for ensemble and typical averaging. We 
used A = 0.125, as this value is compatible with several 
of the system sizes listed in Tab.|lj The fixed A method is 
more stable, since the number of parameters to fit grows 
only linearly with the expansion orders. Stability was a 
serious issue, because the largest system size available, 
L = 56, was only about half of the one used in Refs. [22b 
[M1 Due to this limitation, fits were stable for adding 
or removing an expansion parameter only in the range 
0 < g < 1. The reason for this is that, for large |g|, the 
g-th power in Eq. (H strongly enhances the contribution 
of the few spatial points with very large (if g > 0) or very 
small (if g < 0) wave-function amplitude squared, which 
therefore dominate the sum. This results in an effectively 
reduced statistics and so in a noisy dataset, and leads to 
a regime 0 < g < 1, where GMFEs behave numerically 
the best. Notice that, by construction, Dq = d = 3 and 
Di = ai, and moreover ao .5 = d due to a symmetry 
relation derived in Ref. [33 

The resulting critical parameters are listed in Tab. |TT] 
and shown in Fig. The results are essentially indepen¬ 
dent of g and the type of averaging, as expected. We 
also checked that the critical parameters do not depend 
on the width of the energy window, AE, used in the 
computation of the GMFEs. As we show in Fig. the 
results for E(. and v are independent of AE within errors. 
Moreover, the choice AE = 0.0012 is optimal, as it leads 
to the best accuracy. 

To quote a final result for the critical parameters, we 
have averaged the values of Ac, v and y, and of the corre¬ 
sponding errors, obtained with the various GMFEs. (A 
weighted average, using the inverse of the error band as 
weight, yields similar numbers.) Our result for the crit¬ 
ical point, Ec = 0.3357 (0.3340..0.3368), is compatible 
with the value reported in Ref. |T3 at the 2-a level. On 
average, the systematic error on Ec is = 0.0002, so 


negligible compared to the statistical error. Our result for 
the critical exponent, i/ = 1.461 (1.429..1.519), agrees at 
the 1-cr level with the result of Ref. m and with previous 
result s for the critical exponent of the unitary Anderson 
model2212SI. For this quantity, one has also to take into 
account that on average the systematic error due to trun¬ 
cation, = 0.040, is of the same size as the statistical 
error. On the other hand, our value for the irrelevant ex¬ 
ponent, y = 3.307 (2.210..4.572), is significantly different 
from the value of Ref. 1331 y^^ = 1.651 (1.601..1.707). It 
is well known that it is very difficult to determine irrel¬ 
evant exponents accurately, and to explain this discrep¬ 
ancy further work and higher-quality data are needed. It 
is possible that for the system sizes presently available, 
more than one irrelevant term gives important contribu¬ 
tions, so that our result for y would be a sort of “effec¬ 
tive” irrelevant exponent. In any case this point requires 
further analysis. 

As a final remark, we note that since results for differ¬ 
ent g are strongly correlated, there is no significance in 
the fact that our values for the critical point are system¬ 
atically lower, and the ones for the critical exponent are 
systematically higher than the reference values. 

The convergence of the fixed-A MFSS confirms the 
presence of a critical point in the QGD Dirac spec¬ 
trum where the system undergoes a true localization- 
delocalization transition, employing completely different 
observables than the ones used in Ref. M The results of 
our analysis also provide further evidence that the tran¬ 
sition in the QCD Dirac spectrum belongs to the univer¬ 
sality class of the 3D unitary Anderson model. Moreover, 
despite the fact that it does not provide the values of the 
MFEs, the convergence of this method also strongly indi¬ 
cates the presence of multifractality at the critical point. 

We next procedeed to apply the variable-A method, in 
order to try and determine the multifractal exponents, 
and compare them to the ones obtained for the unitary 


























FIG. 4. Critical parameters obtained via MFFS at fixed A = 0.125 on the eigenvectors of the QCD Dirac operator. Error bars 
correspond to the 95% confidence band. Systematic errors are not included. 
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TABLE II. Result of the MFSS at fixed A = 0.125 for the eigenvectors of the Dirac operator of QCD. 
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FIG. 5. Dependence of the fitted critical point and critical 
exponent, as obtained from Dq'Y’ at fixed A = 0.125, for var¬ 
ious energy windows AE. Error bars correspond to the 95% 
confidence band. Only statistical errors are shown. 


Anderson model. However, this method requires small 
values of A to work properly, and is also more demanding 
as it is a two-variable fit. In practice, the /Ndf ratio 
reached a value close to unity only if we left out the 
smallest system sizes, below Lmin = 36, and if we used 
data corresponding to £ = 1 and 2 only. Although using 
I-'min = 36 and £ = 1,2 improved the convergence, the 


fits were still unstable against changing the expansion 
orders. This can be understood, as a similar amount of 
independent data is available as in the fixed-A method, 
but there are many more parameters to fit, as discussed 
in Sec. |IT] In order to be able to estimate the MFEs, we 
then fixed the critical energy and the critical exponent 
to the values obtained with the fixed-A method, Ec = 
0.3357 and v = 1.461, in this way stabilizing the fits. The 
systematic uncertainty corresponding to this procedure 
was estimated by repeating the fits with and v fixed 
to one of the four possible combinations of the values 
and which are the lower and upper boundaries 
of the confidence interval of Ec and v, respectively. The 
largest and smallest values obtained in this way were then 
used as upper and lower error bar on the MFEs. We 
experienced that the main source of uncertainty comes 
from the choice of Ec, while fits are much less sensitive 
to the choice of v. Moreover, statistical errors (estimated 
by Monte-Carlo) and systematic errors due to truncation 
were comparatively negligible. 

The results of this procedure are depicted in Eig. 

A set of nontrivial MFEs was obtained, thus providing 
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FIG. 6. Estimated values of the MFEs, (a) aq and (b) Dq, in i 
model taken from Ref. |23] (slightly shifted horizontally for clari 
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L-temperature QCD, and MFEs of the 3D unitary Anderson 


direct evidence of the multifractality of the critical eigen¬ 
functions of the QCD Dirac operator. Moreover, our re¬ 
sults for the MFEs in QCD are compatible with the ones 
obtained in the unitary Anderson model, which further 
confirms that the transition belongs to the chiral unitary 
Anderson class. 


VI. SUMMARY 

We investigated the Anderson transition in the spec¬ 
trum of the Dirac operator of QCD at high tempera¬ 
ture, found by the authors of Ref. [TH While that work 
made use of spectral statistics, our aim in this paper 
was to examine the transition by studying the eigenvec¬ 
tors, and their multifractal properties at the critical en¬ 
ergy. The results of Ref. [14] for the correlation length 
critical exponent suggested that the Anderson transition 
in QCD belongs to the same universality class as the 
three-dimensional unitary Anderson model. We there¬ 
fore looked for more similarities between these models. 

First we examined the correlations between eigenvec¬ 
tors of a given gauge configuration. We found strong 
correlations between eigenmodes of the QCD Dirac op¬ 
erator, decreasing with energy separation in a similar way 
as in the unitary Anderson model. We then performed 
two multifractal finite-size scaling (MESS) analyses, one 
with fixed ratio A of the coarse-graining box size to the 
system size, and one with variable A. MESS with the 
fixed-A method allowed an alternative determination of 
the critical point and of the critical exponent, which is 
in agreement with the findings of Ref. [TU and, for the 


critical exponent, with those of Refs. US] and [321 for the 
unitary Anderson model. To perform MESS with the 
variable-A method and determine the multifractal expo¬ 
nents (MFEs), we performed fits fixing the critical energy 
and the critical exponent to the values obtained with the 
fixed-A method. The resulting MFEs are compatible with 
the MFEs found in the unitary Anderson model. 

In conclusion, our work confirms the presence of an An¬ 
derson metal-insulator phase transition in the spectrum 
of the Dirac operator in high-temperature QCD, and pro¬ 
vides further evidence that this transition belongs to the 
three-dimensional unitary Anderson model class. Mor- 
ever, we have shown that the critical wave-functions of 
the Dirac operator are multifractals. The physical con¬ 
sequences of the QCD Anderson transition and of multi¬ 
fractality still largely need to be explored, and may lead 
in particular to a better understanding of the QCD chi¬ 
ral transition. Further work along these lines might prove 
beneficial for condensed matter physics as well, as it ap¬ 
proaches the subject of localization/delocalization tran¬ 
sitions from a broader perspective. 
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